Read in Data

Glasser regions list

glasser.regions <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/Maps/HCPMMP_glasseratlas/glasser360_regionlist.csv")
glasser.snr.exclude <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/Maps/SNR/glasser_SNR_exclusion.csv")
glasser.frontal <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/Maps/HCPMMP_glasseratlas/glasser360_regionlist_frontallobe.csv")
glasser.lh.frontallobe.order <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/Maps/Brainsmash/LeftDistmat_frontallobe_regionorder.txt", header = F) %>% set_names("orig_parcelname")
glasser.rh.frontallobe.order <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/Maps/Brainsmash/RightDistmat_frontallobe_regionorder.txt", header = F) %>% set_names("orig_parcelname")

S-A Axis

SAaxis.cifti <- read_cifti("/Volumes/Hera/Projects/corticalmyelin_development/Maps/S-A_ArchetypalAxis/FSLRVertex/SensorimotorAssociation_Axis_parcellated/SensorimotorAssociation.Axis.Glasser360.pscalar.nii")
SAaxis <- data.frame(SA.axis = rank(SAaxis.cifti$data), orig_parcelname = names(SAaxis.cifti$Parcel))

BigBrain cytoarchitectural variation

bigbrain.cifti <- read_cifti("/Volumes/Hera/Projects/corticalmyelin_development/Maps/BigBrain_histologygradient/BigBrain.Histology.pscalar.nii")
bigbrain <- data.frame(cytoarchitecture.gradient = bigbrain.cifti$data, orig_parcelname = names(bigbrain.cifti$Parcel))

Myelin development data for correlation significance testing

regional.rate <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/code/corticalmyelin_maturation/results/R1_depth_development/regional_averagederivative_frontallobe.csv")
regional.maturationage <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/code/corticalmyelin_maturation/results/R1_depth_development/regional_maturationage_frontallobe.csv")

S-A Axis BrainSmash

Autocorrelation-preserving surrogates

#Hemisphere-specific surrogates generated by brainsmash (n = 500 per hemisphere)
surrogate_maps_rh <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/Maps/Brainsmash/smash_data/SAaxis.rh.frontallobe.surrogates.txt") %>% t() %>% as.data.frame() #72 rows are RH regions, 500 columns are 500 surrogate maps 
surrogate_maps_lh <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/Maps/Brainsmash/smash_data/SAaxis.lh.frontallobe.surrogates.txt") %>% t() %>% as.data.frame() #72 rows are LH regions, 500 columns are 500 surrogate maps 

#Reflect data across hemisphere to generate frontal lobe surrogate maps
surrogate_maps_rh <- rbind(surrogate_maps_rh, surrogate_maps_rh)
surrogate_maps_lh <- rbind(surrogate_maps_lh, surrogate_maps_lh)

#Add in region names in correct order
glasser.frontallobe.order <- rbind(glasser.rh.frontallobe.order, glasser.lh.frontallobe.order)
surrogate_maps_rh$orig_parcelname <- glasser.frontallobe.order$orig_parcelname
surrogate_maps_lh$orig_parcelname <- glasser.frontallobe.order$orig_parcelname

#NA out low SNR parcels to exclude from null correlation tests
surrogate_maps_rh <- surrogate_maps_rh[!(surrogate_maps_rh$orig_parcelname %in% glasser.snr.exclude$orig_parcelname),]
surrogate_maps_rh <- left_join(glasser.frontallobe.order, surrogate_maps_rh, by = "orig_parcelname")

surrogate_maps_lh <- surrogate_maps_lh[!(surrogate_maps_lh$orig_parcelname %in% glasser.snr.exclude$orig_parcelname),]
surrogate_maps_lh <- left_join(glasser.frontallobe.order, surrogate_maps_lh, by = "orig_parcelname")

surrogate_maps <- merge(surrogate_maps_rh, surrogate_maps_lh, by = c("orig_parcelname"), sort = F) #1000 nulls!
colnames(surrogate_maps)[2:1001] <-  sprintf("map%s",seq(from = 1, to = 1000)) 
surrogate_maps <- left_join(surrogate_maps, glasser.regions, by = c("orig_parcelname"))

Visualize Example Brainmash Autocorrelation-Preserving Nulls

for(map in c("map1", "map100", "map200", "map300", "map400", "map500", "map600", "map700", "map800", "map900", "map100")){
  null.map <- surrogate_maps %>% filter(label != "lh_L_10pp") %>% select(label, map)
  colnames(null.map)<- c("label", "nullmap")
  SA.nullmap.plot <- ggseg(.data = null.map, atlas = "glasser", mapping = aes(fill = nullmap), colour=I("gray50"), size=I(.06)) +
  scale_fill_gradient2(low= "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "fill")
  print(SA.nullmap.plot)
}

Compute Brainsmash Null-Based Correlation Significance

#Function to compute the brainsmash-based permutation p-value (pSMASH!) given an empirical correlation between depth-specific development map and the SA-preserving null S-A axis maps #super smash bros
brainsmash.results <- function(dev.df, dev.metric, this.depth){
  df <- dev.df %>% filter(depth == this.depth)
  df <- df %>% select(label, any_of(dev.metric))
  surrogate_maps.dev <- left_join(surrogate_maps, df, by = c("label")) #ensure surrogate map and dev df regional ordering match
  surrogate_maps.dev <- left_join(surrogate_maps.dev, SAaxis, by = c("orig_parcelname"))
  
  surrogate_maps.null <- surrogate_maps.dev %>% select(contains("map")) #1000 null maps
  surrogate_maps.empirical <- surrogate_maps.dev %>% select(SA.axis, any_of(dev.metric)) #true data
  
  empirical.cor <- cor(surrogate_maps.empirical$SA.axis, surrogate_maps.empirical[dev.metric], method = c("spearman"), use = "complete.obs") #true correlation 
  
  null.cors <- cor(surrogate_maps.empirical[dev.metric], surrogate_maps.null, method = c("spearman"), use = "complete.obs") #correlation with brainsmash based surrogates 
  
  mean.nullcor <- mean(null.cors)
  if(empirical.cor > 0){
    brainsmash.p <- sum(null.cors > empirical.cor[1])/1000
  }
  if(empirical.cor < 0){
    brainsmash.p <- sum(null.cors < empirical.cor[1])/1000
  }
  
  smash.results <- list(empirical.cor[1], mean.nullcor, brainsmash.p)
  names(smash.results) <- list("empirical.cor", "average.null.cor", "brainsmash.pvalue")
  return(smash.results)
}

R1 Rate of Increase - S-A Axis Correlation by Depth

brainsmash.results(dev.df = regional.rate, dev.metric = "rate", this.depth = "depth_1")
## $empirical.cor
## [1] -0.5096381
## 
## $average.null.cor
## [1] 0.01730044
## 
## $brainsmash.pvalue
## [1] 0.046
brainsmash.results(dev.df = regional.rate, dev.metric = "rate", this.depth = "depth_7")
## $empirical.cor
## [1] -0.2408838
## 
## $average.null.cor
## [1] 0.02411319
## 
## $brainsmash.pvalue
## [1] 0.254

BigBrain BrainSmash

Autocorrelation-preserving surrogates

#Hemisphere-specific surrogates generated by brainsmash (n = 500 per hemisphere)
surrogate_maps_rh <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/Maps/Brainsmash/smash_data/Bigbrain.rh.frontallobe.surrogates.txt") %>% t() %>% as.data.frame() #72 rows are RH regions, 500 columns are 500 surrogate maps 
surrogate_maps_lh <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/Maps/Brainsmash/smash_data/Bigbrain.lh.frontallobe.surrogates.txt") %>% t() %>% as.data.frame() #72 rows are LH regions, 500 columns are 500 surrogate maps 

#Reflect data across hemisphere to generate frontal lobe surrogate maps
surrogate_maps_rh <- rbind(surrogate_maps_rh, surrogate_maps_rh)
surrogate_maps_lh <- rbind(surrogate_maps_lh, surrogate_maps_lh)

#Add in region names in correct order
glasser.frontallobe.order <- rbind(glasser.rh.frontallobe.order, glasser.lh.frontallobe.order)
surrogate_maps_rh$orig_parcelname <- glasser.frontallobe.order$orig_parcelname
surrogate_maps_lh$orig_parcelname <- glasser.frontallobe.order$orig_parcelname

#NA out low SNR parcels to exclude from null correlation tests
surrogate_maps_rh <- surrogate_maps_rh[!(surrogate_maps_rh$orig_parcelname %in% glasser.snr.exclude$orig_parcelname),]
surrogate_maps_rh <- left_join(glasser.frontallobe.order, surrogate_maps_rh, by = "orig_parcelname")

surrogate_maps_lh <- surrogate_maps_lh[!(surrogate_maps_lh$orig_parcelname %in% glasser.snr.exclude$orig_parcelname),]
surrogate_maps_lh <- left_join(glasser.frontallobe.order, surrogate_maps_lh, by = "orig_parcelname")

surrogate_maps <- merge(surrogate_maps_rh, surrogate_maps_lh, by = c("orig_parcelname"), sort = F) #1000 nulls!
colnames(surrogate_maps)[2:1001] <-  sprintf("map%s",seq(from = 1, to = 1000)) 
surrogate_maps <- left_join(surrogate_maps, glasser.regions, by = c("orig_parcelname"))

Visualize Example Brainmash Autocorrelation-Preserving Nulls

for(map in c("map1", "map100", "map200", "map300", "map400", "map500", "map600", "map700", "map800", "map900", "map100")){
  null.map <- surrogate_maps %>% filter(label != "lh_L_10pp") %>% select(label, map)
  colnames(null.map)<- c("label", "nullmap")
  BB.nullmap.plot <- ggseg(.data = null.map, atlas = "glasser", mapping = aes(fill = nullmap), colour=I("gray50"), size=I(.06)) +
  scale_fill_gradient2(low= "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "fill")
  print(BB.nullmap.plot)
}

Compute Brainsmash Null-Based Correlation Significance

#Function to compute the brainsmash-based permutation p-value (pSMASH!) given an empirical correlation between depth-specific development map and the SA-preserving null BigBrain maps #big brain smash
brainsmash.results <- function(dev.df, dev.metric, this.depth){
  df <- dev.df %>% filter(depth == this.depth)
  df <- df %>% select(label, any_of(dev.metric))
  surrogate_maps.dev <- left_join(surrogate_maps, df, by = c("label")) #ensure surrogate map and dev df regional ordering match
  surrogate_maps.dev <- left_join(surrogate_maps.dev, bigbrain, by = c("orig_parcelname"))
  
  surrogate_maps.null <- surrogate_maps.dev %>% select(contains("map")) #1000 null maps
  surrogate_maps.empirical <- surrogate_maps.dev %>% select(cytoarchitecture.gradient, any_of(dev.metric)) #true data
  
  empirical.cor <- cor(surrogate_maps.empirical$cytoarchitecture.gradient, surrogate_maps.empirical[dev.metric], method = c("spearman"), use = "complete.obs") #true correlation 
  
  null.cors <- cor(surrogate_maps.empirical[dev.metric], surrogate_maps.null, method = c("spearman"), use = "complete.obs") #correlation with brainsmash based surrogates 
  
  mean.nullcor <- mean(null.cors)
  if(empirical.cor > 0){
    brainsmash.p <- sum(null.cors > empirical.cor[1])/1000
  }
  if(empirical.cor < 0){
    brainsmash.p <- sum(null.cors < empirical.cor[1])/1000
  }
  
  smash.results <- list(empirical.cor[1], mean.nullcor, brainsmash.p)
  names(smash.results) <- list("empirical.cor", "average.null.cor", "brainsmash.pvalue")
  return(smash.results)
}

R1 Rate of Increase - BigBrain Cytoarchitectural Variation Correlation by Depth

brainsmash.results(dev.df = regional.rate, dev.metric = "rate", this.depth = "depth_1")
## $empirical.cor
## [1] 0.6371086
## 
## $average.null.cor
## [1] -0.004780919
## 
## $brainsmash.pvalue
## [1] 0.013
brainsmash.results(dev.df = regional.rate, dev.metric = "rate", this.depth = "depth_7")
## $empirical.cor
## [1] 0.5435518
## 
## $average.null.cor
## [1] -0.004175066
## 
## $brainsmash.pvalue
## [1] 0.1